In this tutorial we compare the formulation of a Navier-Stokes problem using standard assembly with mixed function spaces and block assembly.
import typing
import dolfinx.fem
import dolfinx.fem.petsc
import dolfinx.io
import gmsh
import mpi4py.MPI
import numpy as np
import numpy.typing
import petsc4py.PETSc
import ufl
import viskex
import multiphenicsx.fem
import multiphenicsx.fem.petsc
nu = 0.01
def u_in_eval(x: np.typing.NDArray[np.float64]) -> np.typing.NDArray[ # type: ignore[no-any-unimported]
petsc4py.PETSc.ScalarType]:
"""Return the flat velocity profile at the inlet."""
values = np.zeros((2, x.shape[1]))
values[0, :] = 1.0
return values
def u_wall_eval(x: np.typing.NDArray[np.float64]) -> np.typing.NDArray[ # type: ignore[no-any-unimported]
petsc4py.PETSc.ScalarType]:
"""Return the zero velocity at the wall."""
return np.zeros((2, x.shape[1]))
pre_step_length = 4.
after_step_length = 14.
pre_step_height = 3.
after_step_height = 5.
lcar = 1. / 5.
gmsh.initialize()
gmsh.model.add("mesh")
p0 = gmsh.model.geo.addPoint(0.0, after_step_height - pre_step_height, 0.0, lcar)
p1 = gmsh.model.geo.addPoint(pre_step_length, after_step_height - pre_step_height, 0.0, lcar)
p2 = gmsh.model.geo.addPoint(pre_step_length, 0.0, 0.0, lcar)
p3 = gmsh.model.geo.addPoint(pre_step_length + after_step_length, 0.0, 0.0, lcar)
p4 = gmsh.model.geo.addPoint(pre_step_length + after_step_length, after_step_height, 0.0, lcar)
p5 = gmsh.model.geo.addPoint(0.0, after_step_height, 0.0, lcar)
l0 = gmsh.model.geo.addLine(p0, p1)
l1 = gmsh.model.geo.addLine(p1, p2)
l2 = gmsh.model.geo.addLine(p2, p3)
l3 = gmsh.model.geo.addLine(p3, p4)
l4 = gmsh.model.geo.addLine(p4, p5)
l5 = gmsh.model.geo.addLine(p5, p0)
line_loop = gmsh.model.geo.addCurveLoop([l0, l1, l2, l3, l4, l5])
domain = gmsh.model.geo.addPlaneSurface([line_loop])
gmsh.model.geo.synchronize()
gmsh.model.addPhysicalGroup(1, [l5], 1)
gmsh.model.addPhysicalGroup(1, [l0, l1, l2, l4], 2)
gmsh.model.addPhysicalGroup(2, [domain], 0)
gmsh.model.mesh.generate(2)
Info : Meshing 1D... Info : [ 0%] Meshing curve 1 (Line) Info : [ 20%] Meshing curve 2 (Line) Info : [ 40%] Meshing curve 3 (Line) Info : [ 50%] Meshing curve 4 (Line) Info : [ 70%] Meshing curve 5 (Line) Info : [ 90%] Meshing curve 6 (Line) Info : Done meshing 1D (Wall 0.000660553s, CPU 0.000767s) Info : Meshing 2D... Info : Meshing surface 1 (Plane, Frontal-Delaunay) Info : Done meshing 2D (Wall 0.0681941s, CPU 0.068262s) Info : 2520 nodes 5044 elements
mesh, subdomains, boundaries = dolfinx.io.gmshio.model_to_mesh(
gmsh.model, comm=mpi4py.MPI.COMM_WORLD, rank=0, gdim=2)
gmsh.finalize()
boundaries_1 = boundaries.indices[boundaries.values == 1]
boundaries_2 = boundaries.indices[boundaries.values == 2]
viskex.dolfinx.plot_mesh(mesh)
viskex.dolfinx.plot_mesh_tags(mesh, boundaries, "boundaries")
viskex.dolfinx.plot_mesh_entities(mesh, mesh.topology.dim - 1, "boundaries_1", boundaries_1)
viskex.dolfinx.plot_mesh_entities(mesh, mesh.topology.dim - 1, "boundaries_2", boundaries_2)
V_element = ufl.VectorElement("Lagrange", mesh.ufl_cell(), 2)
Q_element = ufl.FiniteElement("Lagrange", mesh.ufl_cell(), 1)
def run_monolithic() -> dolfinx.fem.Function:
"""Run standard FEniCSx formulation using a mixed function space."""
# Function spaces
W_element = ufl.MixedElement(V_element, Q_element)
W = dolfinx.fem.FunctionSpace(mesh, W_element)
# Test and trial functions: monolithic
vq = ufl.TestFunction(W)
(v, q) = ufl.split(vq)
dup = ufl.TrialFunction(W)
up = dolfinx.fem.Function(W)
(u, p) = ufl.split(up)
# Variational forms
F = (nu * ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx
+ ufl.inner(ufl.grad(u) * u, v) * ufl.dx
- ufl.inner(p, ufl.div(v)) * ufl.dx
+ ufl.inner(ufl.div(u), q) * ufl.dx)
J = ufl.derivative(F, up, dup)
# Boundary conditions
u_in = dolfinx.fem.Function(W.sub(0).collapse()[0])
u_in.interpolate(u_in_eval)
u_wall = dolfinx.fem.Function(W.sub(0).collapse()[0])
u_wall.interpolate(u_wall_eval)
bdofs_V_1 = dolfinx.fem.locate_dofs_topological(
(W.sub(0), W.sub(0).collapse()[0]), mesh.topology.dim - 1, boundaries_1)
bdofs_V_2 = dolfinx.fem.locate_dofs_topological(
(W.sub(0), W.sub(0).collapse()[0]), mesh.topology.dim - 1, boundaries_2)
inlet_bc = dolfinx.fem.dirichletbc(u_in, bdofs_V_1, W.sub(0))
wall_bc = dolfinx.fem.dirichletbc(u_wall, bdofs_V_2, W.sub(0))
bc = [inlet_bc, wall_bc]
# Class for interfacing with SNES
class NavierStokesProblem(object):
"""Define a nonlinear problem, interfacing with SNES."""
def __init__( # type: ignore[no-any-unimported]
self, F: ufl.Form, J: ufl.Form, solution: dolfinx.fem.Function,
bcs: typing.List[dolfinx.fem.DirichletBC], P: typing.Optional[ufl.Form] = None
) -> None:
self._F = dolfinx.fem.form(F)
self._J = dolfinx.fem.form(J)
self._obj_vec = dolfinx.fem.petsc.create_vector(self._F)
self._solution = solution
self._bcs = bcs
self._P = P
def create_snes_solution(self) -> petsc4py.PETSc.Vec: # type: ignore[no-any-unimported]
"""
Create a petsc4py.PETSc.Vec to be passed to petsc4py.PETSc.SNES.solve.
The returned vector will be initialized with the initial guess provided in `self._solution`.
"""
x = self._solution.vector.copy()
with x.localForm() as _x, self._solution.vector.localForm() as _solution:
_x[:] = _solution
return x
def update_solution(self, x: petsc4py.PETSc.Vec) -> None: # type: ignore[no-any-unimported]
"""Update `self._solution` with data in `x`."""
x.ghostUpdate(addv=petsc4py.PETSc.InsertMode.INSERT, mode=petsc4py.PETSc.ScatterMode.FORWARD)
with x.localForm() as _x, self._solution.vector.localForm() as _solution:
_solution[:] = _x
def obj( # type: ignore[no-any-unimported]
self, snes: petsc4py.PETSc.SNES, x: petsc4py.PETSc.Vec
) -> np.float64:
"""Compute the norm of the residual."""
self.F(snes, x, self._obj_vec)
return self._obj_vec.norm() # type: ignore[no-any-return]
def F( # type: ignore[no-any-unimported]
self, snes: petsc4py.PETSc.SNES, x: petsc4py.PETSc.Vec, F_vec: petsc4py.PETSc.Vec
) -> None:
"""Assemble the residual."""
self.update_solution(x)
with F_vec.localForm() as F_vec_local:
F_vec_local.set(0.0)
dolfinx.fem.petsc.assemble_vector(F_vec, self._F)
dolfinx.fem.apply_lifting(F_vec, [self._J], [self._bcs], x0=[x], scale=-1.0)
F_vec.ghostUpdate(addv=petsc4py.PETSc.InsertMode.ADD, mode=petsc4py.PETSc.ScatterMode.REVERSE)
dolfinx.fem.set_bc(F_vec, self._bcs, x, -1.0)
def J( # type: ignore[no-any-unimported]
self, snes: petsc4py.PETSc.SNES, x: petsc4py.PETSc.Vec, J_mat: petsc4py.PETSc.Mat,
P_mat: petsc4py.PETSc.Mat
) -> None:
"""Assemble the jacobian."""
J_mat.zeroEntries()
dolfinx.fem.petsc.assemble_matrix( # type: ignore[misc]
J_mat, self._J, self._bcs, diagonal=1.0) # type: ignore[arg-type]
J_mat.assemble()
if self._P is not None:
P_mat.zeroEntries()
dolfinx.fem.petsc.assemble_matrix( # type: ignore[misc]
P_mat, self._P, self._bcs, diagonal=1.0) # type: ignore[arg-type]
P_mat.assemble()
# Create problem
problem = NavierStokesProblem(F, J, up, bc)
F_vec = dolfinx.fem.petsc.create_vector(problem._F)
J_mat = dolfinx.fem.petsc.create_matrix(problem._J)
# Solve
snes = petsc4py.PETSc.SNES().create(mesh.comm)
snes.setTolerances(max_it=20)
snes.getKSP().setType("preonly")
snes.getKSP().getPC().setType("lu")
snes.getKSP().getPC().setFactorSolverType("mumps")
snes.setObjective(problem.obj)
snes.setFunction(problem.F, F_vec)
snes.setJacobian(problem.J, J=J_mat, P=None)
snes.setMonitor(lambda _, it, residual: print(it, residual))
up_copy = problem.create_snes_solution()
snes.solve(None, up_copy)
problem.update_solution(up_copy) # TODO can this be safely removed?
up_copy.destroy()
snes.destroy()
return up
up_m = run_monolithic()
(u_m, p_m) = (up_m.sub(0).collapse(), up_m.sub(1).collapse())
0 5.423530856245577 1 0.12945659318075842 2 0.04875205308048891 3 0.012455097181191853 4 0.010454482183770904 5 0.0008411904377640517 6 4.2733147396805684e-05 7 5.832046606630726e-08 8 2.8939001202974883e-13
viskex.dolfinx.plot_vector_field(u_m, "u")
viskex.dolfinx.plot_vector_field(u_m, "u", glyph_factor=1.0)